Food competition analysis

Load packages

home <- here::here()

library(tidyverse)
library(tidylog)
library(RCurl)
library(sdmTMB)
library(RColorBrewer)
library(devtools)
library(patchwork)
library(ggstats)
library(ggh4x)
library(viridis)
library(sdmTMBextra)
library(ggcorrplot)

# Source map-plot
#source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
source(paste0(home, "/R/functions/map-plot.R"))

Reading layer `StatRec_map_Areas_Full_20170124' from data source 
  `/Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/shapefiles/ICES-StatRec-mapto-ICES-Areas/StatRec_map_Areas_Full_20170124.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 11074 features and 17 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -44 ymin: 36 xmax: 69 ymax: 85
Geodetic CRS:  WGS 84

Read data & prepare data

d <- read_csv(paste0(home, "/data/clean/aggregated_stomach_data.csv")) |> 
  drop_na(group) |> 
  drop_na(oxy)

# Calculate relative prey weights (saduria and benthos)
d <- d |> 
  drop_na(group) |> 
  drop_na(oxy) |> 
  rename(oxygen = oxy) %>%
  mutate(tot_weight = rowSums(select(., ends_with('_tot'))),  
         benthic_weight = amphipoda_tot + bivalvia_tot + gadus_morhua_tot +
           gobiidae_tot + mysidae_tot + non_bio_tot + 
           other_crustacea_tot + other_tot + other_pisces_tot + platichthys_flesus_tot +
           polychaeta_tot + saduria_entomon_tot) |> 
  rename(saduria_weight = saduria_entomon_tot,
         flounder_density = kg_km2_flounder,
         large_cod_density = kg_km2_large_cod,
         small_cod_density = kg_km2_small_cod) |> 
  mutate(tot_rel_weight = tot_weight / (pred_weight_g - tot_weight), 
         benthic_rel_weight = benthic_weight / (pred_weight_g - tot_weight),
         saduria_rel_weight = saduria_weight / (pred_weight_g - tot_weight)) |> 
  dplyr::select(-ends_with("_tot")) |> 
  dplyr::select(-predator_latin_name, date) |> 
  # Add small constant to density because we want to take the log of it
  mutate(large_cod_density = ifelse(large_cod_density == 0,
                                    min(filter(d, kg_km2_large_cod > 0)$kg_km2_large_cod)*0.5,
                                    large_cod_density),
         flounder_density = ifelse(flounder_density == 0,
                                   min(filter(d, kg_km2_flounder > 0)$kg_km2_flounder)*0.5,
                                   flounder_density),
         small_cod_density = ifelse(small_cod_density == 0,
                                    min(filter(d, kg_km2_small_cod > 0)$kg_km2_small_cod)*0.5,
                                    small_cod_density)) |>
  # scale variables
  mutate(fyear = as.factor(year),
         fquarter = as.factor(quarter),
         fhaul_id = as.factor(haul_id),
         depth_sc = as.numeric(scale(depth)),
         oxygen_sc = as.numeric(scale(oxygen)),
         density_saduria_sc = as.numeric(scale(density_saduria)),
         flounder_density_sc = as.numeric(scale(log(flounder_density))),
         large_cod_density_sc = as.numeric(scale(log(large_cod_density))),
         small_cod_density_sc = as.numeric(scale(log(small_cod_density)))) |> 
  # Scale length by group ..
  mutate(pred_length_cm_sc = as.numeric(scale(pred_length_cm)),
         .by = group)


head(d) |> 
  as.data.frame()

         pred_id saduria_weight species pred_weight_g pred_length_cm year
1   2015_4_COD_1              0     Cod           603             41 2015
2 2015_4_COD_101              0     Cod           317             33 2015
3 2015_4_COD_103              0     Cod           122             23 2015
4 2015_4_COD_104              0     Cod            97             21 2015
5 2015_4_COD_106              0     Cod           442             36 2015
6 2015_4_COD_108              0     Cod           596             42 2015
  quarter month ices_rect subdiv  haul_id smhi_serial_no        X        Y
1       4    11      40G4     25 2015_4_2              1 474.8173 6165.344
2       4    11      40G4     25 2015_4_4              3 468.5617 6170.950
3       4    11      40G4     25 2015_4_4              3 468.5617 6170.950
4       4    11      40G4     25 2015_4_4              3 468.5617 6170.950
5       4    11      40G4     25 2015_4_4              3 468.5617 6170.950
6       4    11      40G4     25 2015_4_4              3 468.5617 6170.950
       lat  lon       date small_cod_density large_cod_density
1 55.63333 14.6 2015-11-20          21.17608         1342.3040
2 55.68333 14.5 2015-11-20          15.88554          344.2489
3 55.68333 14.5 2015-11-20          15.88554          344.2489
4 55.68333 14.5 2015-11-20          15.88554          344.2489
5 55.68333 14.5 2015-11-20          15.88554          344.2489
6 55.68333 14.5 2015-11-20          15.88554          344.2489
  n_catch_small_cod n_catch_large_cod flounder_density n_catch_flounder
1            120.55           1870.30         22.01796             38.0
2            101.90            579.65        513.21620           1058.2
3            101.90            579.65        513.21620           1058.2
4            101.90            579.65        513.21620           1058.2
5            101.90            579.65        513.21620           1058.2
6            101.90            579.65        513.21620           1058.2
      group month_year depth density_saduria   oxygen tot_weight benthic_weight
1 large cod    11_2015 58.35        213.1876 4.605775       0.03           0.03
2 large cod    11_2015 49.50       1421.5291 4.729146       0.13           0.13
3 small cod    11_2015 49.50       1421.5291 4.729146       0.45           0.45
4 small cod    11_2015 49.50       1421.5291 4.729146       0.01           0.01
5 large cod    11_2015 49.50       1421.5291 4.729146       0.00           0.00
6 large cod    11_2015 49.50       1421.5291 4.729146       0.00           0.00
  tot_rel_weight benthic_rel_weight saduria_rel_weight fyear fquarter fhaul_id
1   4.975372e-05       4.975372e-05                  0  2015        4 2015_4_2
2   4.102629e-04       4.102629e-04                  0  2015        4 2015_4_4
3   3.702180e-03       3.702180e-03                  0  2015        4 2015_4_4
4   1.031034e-04       1.031034e-04                  0  2015        4 2015_4_4
5   0.000000e+00       0.000000e+00                  0  2015        4 2015_4_4
6   0.000000e+00       0.000000e+00                  0  2015        4 2015_4_4
    depth_sc oxygen_sc density_saduria_sc flounder_density_sc
1  0.3115294 -0.887890         -0.8350312          -0.6030525
2 -0.3472963 -0.799384          0.1740300           1.7808622
3 -0.3472963 -0.799384          0.1740300           1.7808622
4 -0.3472963 -0.799384          0.1740300           1.7808622
5 -0.3472963 -0.799384          0.1740300           1.7808622
6 -0.3472963 -0.799384          0.1740300           1.7808622
  large_cod_density_sc small_cod_density_sc pred_length_cm_sc
1            1.6627540            0.6180239         1.1002898
2            0.9837031            0.5006654        -0.1718711
3            0.9837031            0.5006654         0.9516973
4            0.9837031            0.5006654         0.4414579
5            0.9837031            0.5006654         0.3051892
6            0.9837031            0.5006654         1.2593099

# d |> 
#   filter(group == "large cod") |> 
#   filter(large_cod_density == 0) |> as.data.frame()
# 2021_1_89

t <- d |> filter(tot_rel_weight == 0)
nrow(t)/nrow(d)

[1] 0.254858

Data overview

d |> 
  rename("Relative benthic weight" = "benthic_rel_weight",
         "Relative Saduria weight" = "saduria_rel_weight") |> 
  pivot_longer(c("Relative benthic weight", "Relative Saduria weight")) |>
  mutate(group = str_to_sentence(group)) |> 
  ggplot(aes(value)) + 
  ggh4x::facet_grid2(factor(group, levels = c("Flounder", "Small cod", "Large cod"))~name,
                     scales = "free", independent = "all") + 
  geom_density(color = NA, alpha = 0.8, fill = "grey30") + 
  labs(y = "Density", x = "Value") + 
  NULL

rename: renamed 2 variables (Relative benthic weight, Relative Saduria weight)

pivot_longer: reorganized (Relative benthic weight, Relative Saduria weight) into (name, value) [was 9366x43, now 18732x43]

mutate: changed 18,732 values (100%) of 'group' (0 new NA)

ggsave(paste0(home, "/figures/supp/data_distribution.pdf"), width = 17, height = 17, units = "cm")

# Size distribution by year
d |> 
  distinct(pred_id, .keep_all = TRUE) |> 
  mutate(group = str_to_sentence(group)) |> 
  ggplot(aes(pred_length_cm, color = as.factor(year))) + 
  facet_wrap(~group, scales = "free", ncol = 1) + 
  geom_density(alpha = 0.8, fill = NA) + 
  labs(y = "Density", x = "Predator length (cm)", color = "Year") + 
  scale_color_viridis(discrete = TRUE, option = "mako") +
  NULL

distinct: no rows removed

mutate: changed 9,366 values (100%) of 'group' (0 new NA)

ggsave(paste0(home, "/figures/supp/predator_sizes.pdf"), width = 11, height = 17, units = "cm")

# Sample size per haul
d |> 
  summarise(n = length(unique(pred_id)), .by = haul_id) |> 
  ggplot(aes(n)) +
  geom_histogram()

summarise: now 350 rows and 2 columns, ungrouped

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# How many with 3 or fewer and what are their sample sizes?
# d |> 
#   summarise(n = length(unique(pred_id)), .by = haul_id) |> 
#   mutate(s = ifelse(n < 3, "s", "l")) |> 
#   summarise(n = n(), .by = s)
# # (29/316)*100
# 
# # What is the size range of those 10%?
# d |> 
#   mutate(n = length(unique(pred_id)), .by = haul_id) |> 
#   filter(n > 1 & n < 11) |> 
#   summarise(max = max(pred_length_cm),
#             min = min(pred_length_cm),
#             .by = c(group, haul_id)) |> 
#   mutate(diff = max - min) |> 
#   as.data.frame()
#   #summarise(mean = mean(diff))

# Calculate weights

# The reason to round here is because data come as per hour
d <- d |> 
  mutate(sampled_n = n(), .by = c(haul_id, group)) |> 
  mutate(f_weight = round(n_catch_flounder) / sampled_n,
         sc_weight = round(n_catch_small_cod) / sampled_n,
         lc_weight = round(n_catch_large_cod) / sampled_n) |> 
  mutate(f_weight = f_weight/mean(f_weight),
         sc_weight = sc_weight/mean(sc_weight),
         lc_weight = lc_weight/mean(lc_weight)) |> 
  mutate(f_weight2 = f_weight*2,
         sc_weight2 = sc_weight*2,
         lc_weight2 = lc_weight*2)

mutate: new variable 'sampled_n' (integer) with 43 unique values and 0% NA

mutate: new variable 'f_weight' (double) with 624 unique values and 0% NA

        new variable 'sc_weight' (double) with 465 unique values and 0% NA

        new variable 'lc_weight' (double) with 523 unique values and 0% NA

mutate: changed 9,307 values (99%) of 'f_weight' (0 new NA)

        changed 8,725 values (93%) of 'sc_weight' (0 new NA)

        changed 9,219 values (98%) of 'lc_weight' (0 new NA)

mutate: new variable 'f_weight2' (double) with 624 unique values and 0% NA

        new variable 'sc_weight2' (double) with 465 unique values and 0% NA

        new variable 'lc_weight2' (double) with 523 unique values and 0% NA

Quick explore

Correlation between variables

# Plot correlation between variables
d_cor <- d |>
  dplyr::select("oxygen_sc", "density_saduria_sc", "flounder_density_sc",
                "large_cod_density_sc", "small_cod_density_sc", "depth_sc")

corr <- round(cor(d_cor), 1)

ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
  theme_classic() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.3))

# Sample size 
d |>
  group_by(species) |> 
  summarise(n = n())

group_by: one grouping variable (species)

summarise: now 2 rows and 2 columns, ungrouped

# A tibble: 2 × 2
  species      n
  <chr>    <int>
1 Cod       5484
2 Flounder  3882

d |>
  group_by(species, quarter) |> 
  summarise(n = n())

group_by: 2 grouping variables (species, quarter)

summarise: now 4 rows and 3 columns, one group variable remaining (species)

# A tibble: 4 × 3
# Groups:   species [2]
  species  quarter     n
  <chr>      <dbl> <int>
1 Cod            1  2890
2 Cod            4  2594
3 Flounder       1  2093
4 Flounder       4  1789

Fit models

Groups are: small cod, large cod and flounder. Response variables are: saduria_rel_weight, benthic_rel_weight or total weight. The latter is only for adult cod, because essentially all prey are benthic for small cod and flounder.

Model random effect structure is selected with AIC (see script 02-)

# This is the reason we don't do total weight for flounder and small cod
d |>
  filter(tot_rel_weight > 0) |> 
  group_by(group) |> 
  mutate(ben_prop = benthic_rel_weight / tot_rel_weight) |> 
  summarise(mean_ben_prop = mean(ben_prop))

# A tibble: 3 × 2
  group     mean_ben_prop
  <chr>             <dbl>
1 flounder          0.978
2 large cod         0.614
3 small cod         0.966

Covariates are: ~ 0 + fyear + fquarter + depth_sc + spatial + spatiotemporal random fields + density covariates. For saduria, we use saduria also in interaction with cod and flounder. For cod we use small cod because large and small cod are very correlated. For benthic and total prey, we instead use oxygen, more as a proxy, as the interaction variable

Main models

pred_flounder_sad <- list()
pred_flounder_ben <- list()
pred_cod_sad <- list()
pred_cod_ben <- list()
coef_sad <- list()
coef_ben <- list()
res_sad <- list()
res_ben <- list()
random_sad <- list()
random_ben <- list()
range_sad <- list()
range_ben <- list()

for(i in unique(d$group)) {
  
  dd <- filter(d, group == i)
  
    if(i == "flounder"){
      weigths <- dd$f_weight
      } else if(i == "small cod"){
        weigths <- dd$sc_weight
        } else if(i == "large cod"){
          weigths <- dd$lc_weight
          }
  
  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 5)

  ggplot() +
    inlabru::gg(mesh$mesh) +
    coord_fixed() +
    geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
    labs(x = "Easting (km)", y = "Northing (km)")
  
  ggsave(paste0(home, "/figures/supp/mesh_", i, ".pdf"), width = 17, height = 17, units = "cm")
  
  
  # Saduria model
  
  m_sad <- sdmTMB(saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + pred_length_cm_sc +
                    small_cod_density_sc*density_saduria_sc + 
                    flounder_density_sc*density_saduria_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  weights = weigths, 
                  spatiotemporal = "iid", 
                  spatial = "on",
                  time = "year")
  print(i)
  sanity(m_sad)
  print(m_sad)
  
  
  # Benthic model
  
    if(unique(dd$group) %in% c("large cod", "small cod")) {
      
      
      m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc +
                        small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                      data = dd,
                      mesh = mesh,
                      family = tweedie(),
                      weights = weigths, 
                      spatiotemporal = "iid", 
                      spatial = "off", 
                      time = "year")
      print(i)
      sanity(m_ben)
      print(m_ben)
      
    } else {
      
      m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc +
                        small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                      data = dd,
                      mesh = mesh,
                      family = tweedie(),
                      weights = weigths, 
                      spatiotemporal = "iid", 
                      spatial = "on",
                      time = "year")
      print(i)
      sanity(m_ben)
      print(m_ben)
      
    }
       
   
  # Spatial and spatiotemporal random effects
  d_haul <- dd |>
    distinct(haul_id, .keep_all = TRUE)

  preds_sad <- predict(m_sad, newdata = d_haul)
  preds_ben <- predict(m_ben, newdata = d_haul)

  random_sad[[i]] <- preds_sad
  random_ben[[i]] <- preds_ben

  # Residuals
  samps <- sdmTMBextra::predict_mle_mcmc(m_sad, mcmc_iter = 401, mcmc_warmup = 400)
  mcmc_res <- residuals(m_sad, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_sad[[i]] <- dd

  samps <- sdmTMBextra::predict_mle_mcmc(m_ben, mcmc_iter = 401, mcmc_warmup = 400)
  mcmc_res <- residuals(m_ben, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_ben[[i]] <- dd


  # Ranges
  range_sad[[i]] <- tidy(m_sad, effects = "ran_pars") |> filter(term == "range") |> mutate(group = i, model = "saduria")
  range_ben[[i]] <- tidy(m_ben, effects = "ran_pars") |> filter(term == "range") |> mutate(group = i, model = "benthos")


  # Conditional effects: flounder
  nd_flounder <- data.frame(expand_grid(
    density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
                           mean(d$density_saduria_sc),
                           quantile(d$density_saduria_sc, probs = 0.95)),
    flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
                              quantile(dd$flounder_density_sc, probs = 0.95),
                              length.out = 50))) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           pred_length_cm_sc = 0,
           oxygen_sc = 0,
           depth_sc = 0,
           small_cod_density_sc = 0)

  preds_flounder_sad <- predict(m_sad, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  preds_flounder_ben <- predict(m_ben, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  pred_flounder_sad[[i]] <- preds_flounder_sad |> mutate(group = i, xvar = "flounder")
  pred_flounder_ben[[i]] <- preds_flounder_ben |> mutate(group = i, xvar = "flounder")

  # Conditional effects: cod
  nd_cod <- data.frame(expand_grid(
    density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
                           quantile(d$density_saduria_sc, probs = 0.95)),
    small_cod_density_sc = seq(quantile(dd$small_cod_density_sc, probs = 0.05),
                               quantile(dd$small_cod_density_sc, probs = 0.95),
                               length.out = 50))) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           pred_length_cm_sc = 0,
           oxygen_sc = 0,
           depth_sc = 0,
           flounder_density_sc = 0) #

  preds_cod_sad <- predict(m_sad, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  preds_cod_ben <- predict(m_ben, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  pred_cod_sad[[i]] <- preds_cod_sad |> mutate(group = i, xvar = "cod")
  pred_cod_ben[[i]] <- preds_cod_ben |> mutate(group = i, xvar = "cod")

  # Coefficients
  coefs_sad <- bind_rows(tidy(m_sad, effects = "fixed", conf.int = TRUE)) |>
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

  coefs_ben <- bind_rows(tidy(m_ben, effects = "fixed", conf.int = TRUE)) |>
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

  coef_sad[[i]] <- coefs_sad |> mutate(group = i)
  coef_ben[[i]] <- coefs_ben |> mutate(group = i)

}

filter: removed 5,629 rows (60%), 3,737 rows remaining

Loading required namespace: INLA

[1] "large cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     pred_length_cm_sc + small_cod_density_sc * density_saduria_sc + 
 Formula:     flounder_density_sc * density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                                        coef.est coef.se
fyear2015                                  -8.23    0.91
fyear2016                                 -10.13    0.70
fyear2017                                 -10.02    0.72
fyear2018                                  -8.72    0.80
fyear2019                                 -11.29    1.13
fyear2020                                  -9.85    0.69
fyear2021                                 -10.09    0.74
fyear2022                                 -10.38    0.77
fquarter4                                  -0.37    0.38
depth_sc                                   -0.64    0.35
oxygen_sc                                   0.34    0.28
pred_length_cm_sc                          -0.69    0.09
small_cod_density_sc                        0.21    0.27
density_saduria_sc                          0.59    0.32
flounder_density_sc                        -0.57    0.27
small_cod_density_sc:density_saduria_sc     0.07    0.26
density_saduria_sc:flounder_density_sc      0.23    0.24

Dispersion parameter: 0.08
Tweedie p: 1.43
Matérn range: 17.53
Spatial SD: 2.27
Spatiotemporal IID SD: 1.77
ML criterion at convergence: -671.212

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "large cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     small_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                               coef.est coef.se
fyear2015                         -6.23    0.38
fyear2016                         -7.09    0.25
fyear2017                         -6.47    0.27
fyear2018                         -6.57    0.30
fyear2019                         -7.52    0.41
fyear2020                         -6.51    0.27
fyear2021                         -6.48    0.26
fyear2022                         -6.44    0.28
fquarter4                          0.76    0.13
depth_sc                          -0.04    0.09
pred_length_cm_sc                 -0.25    0.03
small_cod_density_sc              -0.17    0.10
oxygen_sc                          0.49    0.10
flounder_density_sc               -0.04    0.10
small_cod_density_sc:oxygen_sc     0.02    0.06
oxygen_sc:flounder_density_sc     -0.01    0.08

Dispersion parameter: 0.26
Tweedie p: 1.63
Matérn range: 10.08
Spatiotemporal IID SD: 1.47
ML criterion at convergence: -7083.399

See ?tidy.sdmTMB to extract these values as a data frame.

distinct: removed 3,426 rows (92%), 311 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005737 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 57.37 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 74.047 seconds (Warm-up)
Chain 1:                0.215 seconds (Sampling)
Chain 1:                74.262 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005569 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 55.69 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 36.407 seconds (Warm-up)
Chain 1:                0.049 seconds (Sampling)
Chain 1:                36.456 seconds (Total)
Chain 1: 

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

filter: removed 3 rows (75%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

filter: removed 7,619 rows (81%), 1,747 rows remaining

[1] "small cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     pred_length_cm_sc + small_cod_density_sc * density_saduria_sc + 
 Formula:     flounder_density_sc * density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                                        coef.est coef.se
fyear2015                                  -9.15    1.39
fyear2016                                 -10.46    1.17
fyear2017                                 -10.59    1.10
fyear2018                                 -10.74    1.51
fyear2019                                  -9.37    1.48
fyear2020                                 -10.09    1.12
fyear2021                                 -11.39    1.03
fyear2022                                 -11.03    1.14
fquarter4                                  -1.59    0.36
depth_sc                                   -0.95    0.51
oxygen_sc                                  -0.32    0.34
pred_length_cm_sc                           1.13    0.08
small_cod_density_sc                        0.17    0.35
density_saduria_sc                          0.99    0.45
flounder_density_sc                        -1.20    0.34
small_cod_density_sc:density_saduria_sc    -0.57    0.34
density_saduria_sc:flounder_density_sc      0.61    0.34

Dispersion parameter: 0.12
Tweedie p: 1.49
Matérn range: 23.61
Spatial SD: 3.41
Spatiotemporal IID SD: 2.63
ML criterion at convergence: -1701.241

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "small cod"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     small_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                               coef.est coef.se
fyear2015                         -6.02    0.35
fyear2016                         -5.74    0.22
fyear2017                         -5.48    0.21
fyear2018                         -5.47    0.27
fyear2019                         -5.96    0.35
fyear2020                         -5.66    0.22
fyear2021                         -5.99    0.19
fyear2022                         -5.80    0.21
fquarter4                          0.32    0.10
depth_sc                          -0.42    0.09
pred_length_cm_sc                 -0.23    0.02
small_cod_density_sc               0.00    0.09
oxygen_sc                          0.11    0.08
flounder_density_sc               -0.02    0.07
small_cod_density_sc:oxygen_sc    -0.07    0.09
oxygen_sc:flounder_density_sc      0.06    0.08

Dispersion parameter: 0.07
Tweedie p: 1.52
Matérn range: 10.14
Spatiotemporal IID SD: 1.30
ML criterion at convergence: -11124.987

See ?tidy.sdmTMB to extract these values as a data frame.

distinct: removed 1,472 rows (84%), 275 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005812 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 58.12 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 172.234 seconds (Warm-up)
Chain 1:                0.451 seconds (Sampling)
Chain 1:                172.685 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.005337 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 53.37 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 46.469 seconds (Warm-up)
Chain 1:                0.05 seconds (Sampling)
Chain 1:                46.519 seconds (Total)
Chain 1: 

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

filter: removed 3 rows (75%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

filter: removed 5,484 rows (59%), 3,882 rows remaining

[1] "flounder"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc + 
 Formula:     pred_length_cm_sc + small_cod_density_sc * density_saduria_sc + 
 Formula:     flounder_density_sc * density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                                        coef.est coef.se
fyear2015                                  -7.31    1.01
fyear2016                                  -5.41    0.89
fyear2017                                  -8.11    0.85
fyear2018                                  -8.58    0.96
fyear2019                                  -6.83    1.01
fyear2020                                  -7.62    0.86
fyear2021                                  -9.63    1.15
fyear2022                                  -8.45    0.90
fquarter4                                  -1.39    0.30
depth_sc                                   -0.46    0.27
oxygen_sc                                   0.14    0.22
pred_length_cm_sc                           0.15    0.06
small_cod_density_sc                        0.55    0.16
density_saduria_sc                          0.16    0.20
flounder_density_sc                        -0.98    0.19
small_cod_density_sc:density_saduria_sc     0.13    0.14
density_saduria_sc:flounder_density_sc     -0.21    0.19

Dispersion parameter: 0.22
Tweedie p: 1.52
Matérn range: 66.75
Spatial SD: 2.07
Spatiotemporal IID SD: 1.43
ML criterion at convergence: -1086.911

See ?tidy.sdmTMB to extract these values as a data frame.
[1] "flounder"

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

Spatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     small_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                               coef.est coef.se
fyear2015                         -4.73    0.28
fyear2016                         -5.01    0.25
fyear2017                         -4.89    0.19
fyear2018                         -5.35    0.23
fyear2019                         -4.89    0.26
fyear2020                         -4.48    0.17
fyear2021                         -5.33    0.23
fyear2022                         -5.06    0.19
fquarter4                          0.05    0.12
depth_sc                          -0.48    0.11
pred_length_cm_sc                 -0.03    0.02
small_cod_density_sc               0.26    0.07
oxygen_sc                         -0.07    0.08
flounder_density_sc               -0.28    0.08
small_cod_density_sc:oxygen_sc    -0.07    0.05
oxygen_sc:flounder_density_sc      0.07    0.07

Dispersion parameter: 0.14
Tweedie p: 1.49
Matérn range: 15.20
Spatial SD: 0.78
Spatiotemporal IID SD: 0.62
ML criterion at convergence: -4221.382

See ?tidy.sdmTMB to extract these values as a data frame.

distinct: removed 3,616 rows (93%), 266 rows remaining

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.006037 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 60.37 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 282.131 seconds (Warm-up)
Chain 1:                0.954 seconds (Sampling)
Chain 1:                283.085 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.006351 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 63.51 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 47.119 seconds (Warm-up)
Chain 1:                0.24 seconds (Sampling)
Chain 1:                47.359 seconds (Total)
Chain 1: 

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

mutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NA

Now do a separate model for adult cod looking at total prey

  dd <- filter(d, group == "large cod")

filter: removed 5,629 rows (60%), 3,737 rows remaining

  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 5)

  # Total model
  m_tot <- sdmTMB(tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc +
                    large_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  weights = dd$lc_weight,
                  spatiotemporal = "iid", 
                  spatial = "on",
                  time = "year")
  
  sanity(m_tot)

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

  print(m_tot)

Spatiotemporal model fit by ML ['sdmTMB']
Formula: tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + pred_length_cm_sc + 
 Formula:     large_cod_density_sc * oxygen_sc + flounder_density_sc * 
 Formula:     oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
 
Conditional model:
                               coef.est coef.se
fyear2015                         -4.09    0.32
fyear2016                         -4.53    0.22
fyear2017                         -4.22    0.23
fyear2018                         -4.30    0.26
fyear2019                         -5.06    0.36
fyear2020                         -4.66    0.23
fyear2021                         -4.43    0.22
fyear2022                         -4.60    0.24
fquarter4                          0.09    0.11
depth_sc                          -0.04    0.08
pred_length_cm_sc                  0.28    0.03
large_cod_density_sc              -0.37    0.11
oxygen_sc                          0.03    0.09
flounder_density_sc                0.25    0.08
large_cod_density_sc:oxygen_sc     0.04    0.09
oxygen_sc:flounder_density_sc      0.09    0.07

Dispersion parameter: 0.60
Tweedie p: 1.72
Matérn range: 19.04
Spatial SD: 0.19
Spatiotemporal IID SD: 0.76
ML criterion at convergence: -7543.688

See ?tidy.sdmTMB to extract these values as a data frame.

  # Residuals
  samps <- sdmTMBextra::predict_mle_mcmc(m_tot, mcmc_iter = 401, mcmc_warmup = 400)

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.006094 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 60.94 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 401 [  0%]  (Warmup)
Chain 1: Iteration:  40 / 401 [  9%]  (Warmup)
Chain 1: Iteration:  80 / 401 [ 19%]  (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%]  (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%]  (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%]  (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%]  (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%]  (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%]  (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%]  (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%]  (Warmup)
Chain 1: Iteration: 401 / 401 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 59.882 seconds (Warm-up)
Chain 1:                0.113 seconds (Sampling)
Chain 1:                59.995 seconds (Total)
Chain 1: 

  mcmc_res <- residuals(m_tot, type = "mle-mcmc", mcmc_samples = samps)
  dd$res <- as.vector(mcmc_res)

  res_tot <- dd
  
  # Range
  range_tot <- tidy(m_tot, effects = "ran_pars") |> filter(term == "range") |> mutate(group = "large cod", model = "total")

filter: removed 4 rows (80%), one row remaining

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'model' (character) with one unique value and 0% NA

  # Spatial and spatiotemporal random effects
  d_haul <- dd |> 
    distinct(haul_id, .keep_all = TRUE)

distinct: removed 3,426 rows (92%), 311 rows remaining

  preds_tot <- predict(m_tot, newdata = d_haul)
  
  # Coefficients
  coefs_tot <- bind_rows(tidy(m_tot, effects = "fixed", conf.int = TRUE)) |> 
    mutate(species = "Cod (m)",
           response = "Saduria",
           sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
           sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))

mutate: new variable 'species' (character) with one unique value and 0% NA

        new variable 'response' (character) with one unique value and 0% NA

        new variable 'sig' (character) with 2 unique values and 0% NA

  coefs_tot <- coefs_tot |> mutate(group = "large cod")

mutate: new variable 'group' (character) with one unique value and 0% NA

  # Conditional effects: flounder
  nd_flounder <- data.frame(expand_grid(
    flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
                              quantile(dd$flounder_density_sc, probs = 0.95),
                              length.out = 50))) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           pred_length_cm_sc = 0,
           oxygen_sc = 0,
           depth_sc = 0,
           large_cod_density_sc = 0)

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'pred_length_cm_sc' (double) with one unique value and 0% NA

        new variable 'oxygen_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'large_cod_density_sc' (double) with one unique value and 0% NA

  preds_flounder_tot <- predict(m_tot, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
  pred_flounder_tot <- preds_flounder_tot |> mutate(group = "large cod", xvar = "flounder")

mutate: new variable 'group' (character) with one unique value and 0% NA

        new variable 'xvar' (character) with one unique value and 0% NA

Make dataframes

coef_df <- bind_rows(bind_rows(coef_sad) |> mutate(model = "Saduria"),
                     bind_rows(coef_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

coef_df <- coef_df |> bind_rows(coefs_tot |> mutate(model = "Total"))

mutate: new variable 'model' (character) with one unique value and 0% NA

pred_cod_df <- bind_rows(bind_rows(pred_cod_sad) |> mutate(model = "Saduria"),
                         bind_rows(pred_cod_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

pred_flounder_df <- bind_rows(bind_rows(pred_flounder_sad) |> mutate(model = "Saduria"),
                              bind_rows(pred_flounder_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

res_df <- bind_rows(bind_rows(res_sad) |> mutate(model = "Saduria"),
                    bind_rows(res_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

res_df <- res_df |> bind_rows(res_tot |> mutate(model = "Total"))

mutate: new variable 'model' (character) with one unique value and 0% NA

random_df <- bind_rows(bind_rows(random_sad) |> mutate(model = "Saduria"),
                       bind_rows(random_ben) |> mutate(model = "Benthos"))

mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NA

random_df <- random_df |> bind_rows(preds_tot |> mutate(model = "Total"))

mutate: new variable 'model' (character) with one unique value and 0% NA

range_df <- bind_rows(range_tot, bind_rows(range_ben), bind_rows(range_sad))

Plot spatial random effects

random_df <- random_df |>
  mutate(group = str_to_sentence(group))

mutate: changed 2,015 values (100%) of 'group' (0 new NA)

# Saduria
plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria"),
             aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~factor(group, levels = c("Flounder", "Small cod", "Large cod"))) +
  labs(color = "Spatial\nrandom effect") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.width = unit(0.6, "cm"),
        legend.key.height = unit(0.2, "cm"))

filter: removed 1,163 rows (58%), 852 rows remaining

ggsave(paste0(home, "/figures/supp/omega_sad.pdf"), width = 17, height = 9, units = "cm")


# Now do benthos (only for flounder)
plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Flounder"),
             aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ group) +
  labs(color = "Spatial\nrandom effect") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "right",
        legend.direction = "vertical",
        legend.key.width = unit(0.4, "cm"),
        legend.key.height = unit(0.4, "cm"))

filter: removed 1,749 rows (87%), 266 rows remaining

ggsave(paste0(home, "/figures/supp/omega_ben.pdf"), width = 11, height = 11, units = "cm")

Plot spatiotemporal random effects

# Saduria
sad_eps_sc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,740 rows (86%), 275 rows remaining

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

sad_eps_sc

ggsave(paste0(home, "/figures/supp/epsilon_sad_small_cod.pdf"), width = 17, height = 17, units = "cm")

sad_eps_lc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,704 rows (85%), 311 rows remaining

sad_eps_lc

ggsave(paste0(home, "/figures/supp/epsilon_sad_large_cod.pdf"), width = 17, height = 17, units = "cm")

sad_eps_f <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Saduria" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,749 rows (87%), 266 rows remaining

sad_eps_f

ggsave(paste0(home, "/figures/supp/epsilon_sad_flounder.pdf"), width = 17, height = 17, units = "cm")


# Benthos
ben_eps_sc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,740 rows (86%), 275 rows remaining

ben_eps_sc

ggsave(paste0(home, "/figures/supp/epsilon_ben_small_cod.pdf"), width = 17, height = 17, units = "cm")

ben_eps_lc <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,704 rows (85%), 311 rows remaining

ben_eps_lc

ggsave(paste0(home, "/figures/supp/epsilon_ben_large_cod.pdf"), width = 17, height = 17, units = "cm")

ben_eps_f <- plot_map_fc +
  geom_point(data = random_df |> filter(model == "Benthos" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
  scale_color_gradient2() +
  facet_wrap(~ year) +
  labs(color = "Spatiotemporal\nrandom effect") +
  theme(legend.position = c(0.84, 0.16),
        axis.text.x = element_text(angle = 90))

filter: removed 1,749 rows (87%), 266 rows remaining

ben_eps_f

ggsave(paste0(home, "/figures/supp/epsilon_ben_flounder.pdf"), width = 17, height = 17, units = "cm")


# Total
tot_eps <- plot_map_fc +
    geom_point(data = preds_tot, aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
    scale_color_gradient2() +
    facet_wrap(~ year) +
    labs(color = "Spatiotemporal\nrandom effect") +
    theme(legend.position = c(0.84, 0.16),
          axis.text.x = element_text(angle = 90))

tot_eps

Plot range

pal <- (brewer.pal(n = 11, name = "RdYlBu")[c(11, 4, 1)])

range_df |> arrange(estimate)

# A tibble: 7 × 7
  term  estimate std.error conf.low conf.high group     model  
  <chr>    <dbl>     <dbl>    <dbl>     <dbl> <chr>     <chr>  
1 range     10.1      3.00     5.63      18.1 large cod benthos
2 range     10.1      2.38     6.40      16.1 small cod benthos
3 range     15.2      4.16     8.89      26.0 flounder  benthos
4 range     17.5      6.14     8.83      34.8 large cod saduria
5 range     19.0      5.03    11.3       31.9 large cod total  
6 range     23.6      8.66    11.5       48.4 small cod saduria
7 range     66.8     15.9     41.9      106.  flounder  saduria

range_df |> 
  mutate(group = str_to_sentence(group),
         model2 = ifelse(model == "benthos", "Benthic prey", model),
         model2 = ifelse(model == "saduria", "Saduria", model2),
         model2 = ifelse(model == "total", "Total prey", model2)) |> 
  ggplot(aes(model2, estimate, color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) + 
  geom_point(size = 2) +
  geom_hline(yintercept = 5, linetype = 2, alpha = 0.5) +
  scale_color_manual(values = pal) + 
  labs(x = "", y = "Range (km)", color = "Predator") + 
  theme(aspect.ratio = 1,
        legend.position = c(0.86, 0.86)) 

mutate: changed 7 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

ggsave(paste0(home, "/figures/supp/ranges.pdf"), width = 11, height = 11, units = "cm")

Plot residuals

# Plot residuals
res_df |>
  mutate(group = str_to_title(group)) |>
  ggplot(aes(sample = res)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_grid(model ~ group) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)

mutate: changed 22,469 values (100%) of 'group' (0 new NA)

ggsave(paste0(home, "/figures/supp/qq_relative_prey_weight.pdf"), width = 17, height = 17, units = "cm")

Plot coefficients

coef_df$term |> unique()

 [1] "fyear2015"                              
 [2] "fyear2016"                              
 [3] "fyear2017"                              
 [4] "fyear2018"                              
 [5] "fyear2019"                              
 [6] "fyear2020"                              
 [7] "fyear2021"                              
 [8] "fyear2022"                              
 [9] "fquarter4"                              
[10] "depth_sc"                               
[11] "oxygen_sc"                              
[12] "pred_length_cm_sc"                      
[13] "small_cod_density_sc"                   
[14] "density_saduria_sc"                     
[15] "flounder_density_sc"                    
[16] "small_cod_density_sc:density_saduria_sc"
[17] "density_saduria_sc:flounder_density_sc" 
[18] "small_cod_density_sc:oxygen_sc"         
[19] "oxygen_sc:flounder_density_sc"          
[20] "large_cod_density_sc"                   
[21] "large_cod_density_sc:oxygen_sc"         

# Fix some names
coef_df2 <- coef_df |>
  filter(!grepl('year', term)) |>
  filter(!grepl('quarter', term)) |>
  mutate(term = str_remove_all(term, "_sc"),
         term = str_remove_all(term, "density"),
         term = str_replace_all(term, "_", ""),
         term = str_replace_all(term, "geco", "ge co"),
         term = str_replace_all(term, "llco", "ll co"),
         term = str_replace(term, ":", " × "),
         term = str_to_sentence(term),
         group = str_to_sentence(group),
         model2 = ifelse(model == "Saduria", "Prey=Saduria", NA),
         model2 = ifelse(model == "Benthos", "Prey=Benthic prey", model2),
         model2 = ifelse(model == "Total", "Prey=Total prey", model2),
         sig2 = ifelse(sig == "Y", "N", "Y"),
         term = ifelse(term == "Predlengthcm", "Predator length", term),
         group2 = paste0("Predator=", group))

filter: removed 56 rows (49%), 59 rows remaining

filter: removed 7 rows (12%), 52 rows remaining

mutate: changed 52 values (100%) of 'term' (0 new NA)

        changed 52 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

        new variable 'sig2' (character) with 2 unique values and 0% NA

        new variable 'group2' (character) with 3 unique values and 0% NA

p1 <- 
  coef_df2 |>
  filter(model == "Saduria") |> 
  ggplot(aes(estimate, term, color = factor(group, levels = c("Flounder", "Small cod", "Large cod")), alpha = sig2)) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  facet_grid(factor(model2, levels = c("Prey=Saduria", "Prey=Benthic prey")) ~ 
               factor(group2, levels = c("Predator=Flounder", "Predator=Small cod", "Predator=Large cod"))) + 
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal) +
  labs(x = "Estimate", y = "", alpha = "95% CI crossing 0", color = "Model") +
  theme(legend.position = "bottom", 
        legend.direction = "vertical",
        axis.title.x = element_blank(),
        axis.text.y = ggtext::element_markdown()) +
  guides(color = "none",
         alpha = "none") +
  NULL

filter: removed 28 rows (54%), 24 rows remaining

p2 <- 
  coef_df2 |> 
  filter(model == "Benthos") |> 
  ggplot(aes(estimate, term, color = factor(group, levels = c("Flounder", "Small cod", "Large cod")), alpha = sig2)) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  #facet_wrap(~factor(model2, levels = c("Prey=Saduria", "Prey=Benthic prey")), ncol = 1) + 
  facet_grid(factor(model2, levels = c("Prey=Saduria", "Prey=Benthic prey")) ~ 
               factor(group2, levels = c("Predator=Flounder", "Predator=Small cod", "Predator=Large cod"))) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal) +
  labs(x = "Estimate", y = "", alpha = "95% CI crossing 0", color = "Model") +
  theme(legend.position = "bottom", 
        legend.direction = "vertical",
        axis.title.x = element_blank(),
        strip.text.x.top = element_blank(),
        axis.text.y = ggtext::element_markdown()) +
  guides(color = "none",
         alpha = "none") +
  NULL

filter: removed 31 rows (60%), 21 rows remaining

p3 <- 
  coef_df2 |> 
  filter(model == "Total") |> 
  ggplot(aes(estimate, term, color = factor(group, levels = c("Flounder", "Small cod", "Large cod")), alpha = sig2)) +
  geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
  #facet_wrap(~factor(model2), ncol = 1, strip.position = "right") + 
  facet_grid(model2 ~ factor(group2, levels = c("Predator=Flounder", "Predator=Small cod", "Predator=Large cod"))) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 2) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal[3]) +
  labs(x = "Estimate", y = "", alpha = "95% CI crossing 0", color = "Model") +
  theme(legend.position = "left", 
        legend.direction = "vertical",
        strip.text.x.top = element_blank(),
        axis.text.y = ggtext::element_markdown()) +
  guides(color = "none", 
         alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
  NULL

filter: removed 45 rows (87%), 7 rows remaining

p3b <- (plot_spacer() | p3) & plot_layout(widths = c(1, 5.3))

(p1 / p2 / (p3b)) +
  plot_annotation(tag_levels = "a") +
  plot_layout(axes = "collect") & coord_cartesian(xlim = c(min(coef_df2$conf.low), max(coef_df2$conf.high)))

ggsave(paste0(home, "/figures/coefs.pdf"), width = 17, height = 15, units = "cm")

Plot year and quarter coefficients

# Fix some names
coef_df3 <- coef_df |>
  filter(grepl('year', term)) |>
  mutate(term = str_remove_all(term, "fyear"),
         group = str_to_sentence(group),
         term = as.numeric(term),
         model2 = ifelse(model == "Benthos", "Benthic prey", model),
         model2 = ifelse(model == "Saduria", "Saduria", model2),
         model2 = ifelse(model == "Total", "Total prey", model2))

filter: removed 59 rows (51%), 56 rows remaining

mutate: converted 'term' from character to double (0 new NA)

        changed 56 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

ggplot(coef_df3, aes(term, exp(estimate), color = factor(group, levels = c("Flounder", "Small cod", "Large cod")), 
                     fill = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
  facet_wrap(~model2, scales = "free", ncol = 1) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.4, width = 0,
                position = position_dodge(width = 0.4)) +
  scale_color_manual(values = pal) +
  scale_fill_manual(values = pal) +
  labs(x = "Year", y = "Standardized coefficient", color = "") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5, ncol = 3),
         fill = "none") +
  theme(legend.position = c(0.5, 0.99),
        legend.direction = "vertical",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0),
        strip.text.x.top = element_text(angle = 0, hjust = 0)) +
  NULL

ggsave(paste0(home, "/figures/supp/coefs_year.pdf"), width = 11, height = 21, units = "cm")

# Now do quarter
coef_df5 <- coef_df |>
  filter(term %in% c("fquarter4")) |> 
  mutate(group = str_to_sentence(group),
         model2 = ifelse(model == "Benthos", "Benthic prey", model),
         model2 = ifelse(model == "Saduria", "Saduria", model2),
         model2 = ifelse(model == "Total", "Total prey", model2)) |> 
  mutate(sig2 = ifelse(sig == "Y", "N", "Y"))

filter: removed 108 rows (94%), 7 rows remaining

mutate: changed 7 values (100%) of 'group' (0 new NA)

        new variable 'model2' (character) with 3 unique values and 0% NA

mutate: new variable 'sig2' (character) with 2 unique values and 0% NA

ggplot(coef_df5, aes(estimate, model2,
                     alpha = sig2,
                     color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
  geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
                position = position_dodge(width = 0.5)) +
  scale_alpha_manual(values = c(1, 0.4)) +
  scale_color_manual(values = pal) +
  labs(x = "", y = "Quarter 4 effect", alpha = "95% CI crossing 0", color = "Predator") +
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
         alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
  geom_stripped_rows(aes(y = model2), inherit.aes = FALSE) +
  theme(legend.position = "bottom",
        legend.direction = "horizontal",
        legend.box = "horizontal",
        legend.box.spacing = unit(-3, "pt"),
        legend.margin = margin(0, 0, 0, 0)) + 
  coord_cartesian(expand = 0)

ggsave(paste0(home, "/figures/supp/coefs_quarter.pdf"), width = 17, height = 11, units = "cm")

Conditional effects

# Which CI?
# https://www.calculator.net/confidence-interval-calculator.html
pred_df <- bind_rows(pred_cod_df, pred_flounder_df) |>
  mutate(group = str_to_sentence(group),
         sad = ifelse(density_saduria_sc == min(density_saduria_sc), "Low", NA),
         sad = ifelse(density_saduria_sc == max(density_saduria_sc), "High", sad)) |> 
  drop_na(sad)

mutate: changed 1,500 values (100%) of 'group' (0 new NA)

        new variable 'sad' (character) with 3 unique values and 20% NA

drop_na: removed 300 rows (20%), 1,200 rows remaining

pred_df2 <- bind_rows(pred_flounder_df,
                      pred_flounder_tot |> mutate(model = "Total")) |> 
  mutate(group = str_to_sentence(group),
         density_saduria_sc = replace_na(density_saduria_sc,
                                         median(density_saduria_sc, na.rm = TRUE)))

mutate: new variable 'model' (character) with one unique value and 0% NA

mutate: changed 50 values (5%) of 'density_saduria_sc' (50 fewer NA)

        changed 950 values (100%) of 'group' (0 new NA)

# 75% CI!!
ggplot(pred_df |> filter(model == "Saduria" & xvar == "flounder"),
       aes(flounder_density_sc, exp(est), color = sad, fill = sad)) +
  geom_ribbon(aes(ymin = exp(est - 1.150*est_se), ymax = exp(est + 1.150*est_se)),
              alpha = 0.3, color = NA) +
  geom_line() +
  facet_wrap(~factor(group, levels = c("Flounder", "Small cod", "Large cod")),
             scales = "free", 
             ncol = 3) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  labs(x = "Flounder density", y = "Relative saduria weight",
       color = "Saduria", fill = "Saduria") + 
  theme(legend.position = c(0.95, 0.84),
        strip.text.x.top = element_text(angle = 0, hjust = 0)) + 
  scale_x_continuous(breaks = c(-1, 0, 1)) +
  NULL

filter: removed 900 rows (75%), 300 rows remaining

ggsave(paste0(home, "/figures/conditional_saduria_flounder.pdf"), width = 19, height = 7, units = "cm")

Showing conditional effects of oxygen on small cod feeding on benthos

  dd <- filter(d, group == "large cod")

filter: removed 5,629 rows (60%), 3,737 rows remaining

  mesh <- make_mesh(dd,
                    xy_cols = c("X", "Y"),
                    cutoff = 5)

  # Benthic model
  m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + 
                    small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(),
                  spatiotemporal = "iid",
                  spatial = "off",
                  time = "year")

  sanity(m_ben)

✔ Non-linear minimizer suggests successful convergence

✔ Hessian matrix is positive definite

✔ No extreme or very small eigenvalues detected

✔ No gradients with respect to fixed effects are >= 0.001

✔ No fixed-effect standard errors are NA

✔ No standard errors look unreasonably large

✔ No sigma parameters are < 0.01

✔ No sigma parameters are > 100

✔ Range parameter doesn't look unreasonably large

  nd <- data.frame(oxygen = seq(quantile(d$oxygen, probs = 0.05), quantile(d$oxygen, probs = 0.95),
                                length.out = 50)) |>
    mutate(year = 2020,
           fyear = as.factor(2020),
           fquarter = as.factor(1),
           density_saduria_sc = 0,
           flounder_density_sc = 0,
           depth_sc = 0,
           small_cod_density_sc = 0) |>
    mutate(oxygen_sc = (oxygen - mean(d$oxygen)) / sd(d$oxygen))

mutate: new variable 'year' (double) with one unique value and 0% NA

        new variable 'fyear' (factor) with one unique value and 0% NA

        new variable 'fquarter' (factor) with one unique value and 0% NA

        new variable 'density_saduria_sc' (double) with one unique value and 0% NA

        new variable 'flounder_density_sc' (double) with one unique value and 0% NA

        new variable 'depth_sc' (double) with one unique value and 0% NA

        new variable 'small_cod_density_sc' (double) with one unique value and 0% NA

mutate: new variable 'oxygen_sc' (double) with 50 unique values and 0% NA

  p <- predict(m_ben, newdata = nd, re_form = NA, re_form_iid = NA, se_fit = TRUE)

  ggplot(p, aes(oxygen, exp(est))) +
    geom_line() +
    theme_sleek(base_size = 14) +
    geom_hline(yintercept = 0.0011, col = "red") +
    geom_hline(yintercept = 0.0016, col = "red") +
    geom_vline(xintercept = 4.8, col = "red") +
    geom_vline(xintercept = 7.6, col = "red") +
    NULL

((0.0016 - 0.0011) / 0.0011) * 100

[1] 45.45455